Getting started

Adapted from: SOSTOOLS' SOSDEMO1 (See Section 4.1 of SOSTOOLS User's Manual) and Example 2.4 of [PJ08]

P. Parrilo and A. Jadbabaie Approximation of the joint spectral radius using sum of squares. Linear Algebra and its Applications, Elsevier (2008), 428, 2385-2402

using DynamicPolynomials
@polyvar x y
p = 2*x^4 + 2*x^3*y - x^2*y^2 + 5*y^4

\[ 2x^{4} + 2x^{3}y - x^{2}y^{2} + 5y^{4} \]

We need to pick an SDP solver, see here for a list of the available choices. We use SOSModel instead of Model to be able to use the >= syntax for Sum-of-Squares constraints.

using SumOfSquares
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
con_ref = @constraint(model, p >= 0)
optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We see above that the solver found a feasible solution. We now inspect this solution:

q = gram_matrix(con_ref)

GramMatrix{Float64, MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}, Float64, SymMatrix{Float64}}([2.0 1.0 -2.9518518518518517; 1.0 4.903703703703703 4.4408920985006264e-17; -2.9518518518518517 4.4408920985006264e-17 5.0], MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[x², xy, y²]))

We can get the SOS decomposition from the gram matrix as follows:

sosdec = SOSDecomposition(q)

(-1.4043287476933917*x^2 - 0.685602733872361*x*y + 2.123635092514818*y^2)^2 + (-0.02034325143230049*x^2 - 2.1053486943394346*x*y - 0.6931517747987983*y^2)^2 + (0.16567112157245217*x^2 - 0.03405099490006233*x*y + 0.098562725879792*y^2)^2

We now seek for the SOS decomposition of the following polynomial:

p = 4*x^4*y^6 + x^2 - x*y^2 + y^2

\[ 4x^{4}y^{6} - xy^{2} + x^{2} + y^{2} \]

We build the same model as previously with this new polynomial. Here we can use Model instead of SOSModel as we explicitly constrain p to belong to the SOS cone with p in SOSCone().

model = Model(solver)
con_ref = @constraint(model, p in SOSCone())
optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We can query the SOS decomposition directly from the constraint reference as follows:

sos_decomposition(con_ref)

(-1.9748040632461135*x^2*y^3 - 3.0030174111027177e-17*x*y^2 - 0.19297758710767499*x*y + 8.926562404670611e-18*x + 0.8633544111907667*y)^2 + (-5.807692927006204e-16*x^2*y^3 - 1.7452642425911875*x*y^2 + 3.795209026381727e-15*x*y + 0.7962855107812753*x - 5.291246828605798e-16*y)^2 + (0.2579336224358901*x^2*y^3 - 2.9986579153469545e-15*x*y^2 - 1.5486392128975748*x*y + 1.239155439866639e-15*x + 0.24383463418983595*y)^2 + (1.586584599585843e-16*x^2*y^3 + 0.27599821010522463*x*y^2 + 4.2856583217874536e-17*x*y + 0.6049209744419546*x + 3.729397617841436e-16*y)^2 + (-0.18335527863613849*x^2*y^3 + 8.138056831012051e-17*x*y^2 - 0.10009637589813157*x*y + 2.7171760564317086e-16*x - 0.4417735074073052*y)^2

This page was generated using Literate.jl.